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Numerical  Modeling  of  Heterogeneous  High  Explosives 


Christopher  M.  Engelhardt1 

Air  Force  Research  Laboratory,  Eglin  AFB,  FL  32542,  USA 


This  report  contains  a  numerical  algorithm  for  modeling  the  detonation  and  explosion  of 
a  heterogeneous  mixture  of  high  explosive  and  small  metal  particles.  The  simulation 
examines  a  spherical  explosive  design  with  a  mixture  of  nitromethane  as  the  high  explosive 
and  steel  as  the  metal  particles.  The  algorithm  provides  a  computational  model  of  the 
detonation  and  explosion  by  producing  position,  velocity,  and  temperature  profiles  for  the 
metal  particles  over  time.  For  the  gas  phase,  the  algorithm  produces  position,  velocity, 
temperature,  density,  and  pressure  profiles  over  time.  This  is  accomplished  by  taking  into 
account  the  initial  position  and  velocity  profiles  for  the  metal  particles,  a  corresponding 
particle  drag  law,  appropriate  explosive  energy  and  detonation  pressure  inputs,  and  a  blast 
wave  solution  that  governs  the  thermodynamic  state  of  the  gas  phase.  The  behavior  of  the 
solid  particles  and  gas  phase  throughout  the  explosion  is  simulated  by  a  coupled,  two-phase 
algorithm.  The  results  of  the  model  are  compared  against  experimental  data  and  critiqued 
on  a  theoretical  level  as  well.  Recommendations  and  plans  for  improvements  to  the 
algorithm  are  discussed.  This  model  is  intended  to  provide  a  sound  representation  of  the 
detonation  as  well  as  insight  into  the  behavior  of  a  heterogeneous  explosive. 


Nomenclature 


Vi 

= 

detonation  velocity 

Eb 

= 

blast  wave  energy 

Eg 

= 

Gurney  energy 

r 

= 

ratio  of  specific  heats 

Mg,  Mp 

= 

total  mass  of  explosive  and  metal 

a 

= 

speed  of  sound 

R 

= 

radius 

t 

= 

time 

Fd 

= 

drag  force  on  metal  particle 

T) 

= 

function  of  radius 

d 

= 

particle  diameter 

* 

= 

velocity  function 

Q 

= 

drag  coefficient 

f 

= 

pressure  function 

p 

= 

gas  density 

¥ 

= 

density  function 

vg,vp 

= 

gas  and  particle  velocity 

fi 

= 

function  of  specific  heat 

P 

= 

dynamic  viscosity 

p 

= 

pressure 

ag,ap 

= 

gas  and  particle  volume  concentration 

M 

= 

mach  number 

q 

= 

heat  transfer  rate 

s 

= 

function  of  specific  heat 

A 

thermal  conductivity  coefficient 

X 

= 

position 

Cp 

= 

specific  heat 

at 

= 

initial  speed  of  sound 

m 

= 

particle  mass 

w 

= 

function  of  velocity  and  speed  of  sound 

T 

= 

temperature 

Q 

= 

function  of  velocity  and  speed  of  sound 

I.  Introduction 

THERE  have  been  several  theoretical  methods  developed  that  model  the  detonation  and  explosion  of  a  high 
explosive.  The  addition  of  many  small  metal  particles  into  the  high  explosive  to  create  a  heterogeneous  mixture 
has  recently  been  studied.  In  1986,  Baer  and  Nunziato  described  a  two-phase  model  using  reactive  granular 
materials.  Their  study  used  a  continuum  theory  which  considered  compressibility,  entropy  inequality,  and  accounted 
for  mass,  momentum,  and  energy  exchange.1  In  1992,  Lanovets  conducted  a  numeric  study  on  inert  solid  particles 
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contained  in  a  spherical  charge.  This  investigation  considered  a  continuum  of  detonation  gas  combined  with  heavy 
particles.2  Another  study  that  was  performed  by  Patankar  and  Joseph  in  2001  describes  particle  flow  using  an 
Eulerian-Lagrangian  setup  for  the  fluid  phase  and  the  solid  phase  respectively.  The  simulation  by  Patankar  and 
Joseph  emphasized  flexibility  in  forces  being  modeled,  differing  particle  sizes  and  properties,  dense  flow,  and 
conservative  computing  requirements.3  Powers  followed  the  previous  work  and  examined  granular  energetic 
materials  using  a  two-phase  model  in  2004. 4 

Using  numerical  techniques  to  model  the  two-phase  flow  of  a  heterogeneous  high  explosive  requires  a  velocity 
initialization  scheme  to  capture  the  initial  detonation  acceleration  as  well  as  a  blast  wave  solution  for  subsequent 
acceleration.  The  Gurney  Equations  were  derived  in  1943  as  a  simple  way  to  approximate  the  velocity  of  the 
explosive  shell  fragments.  This  was  accomplished  using  an  energy  balance  between  the  kinetic  energy  of  the 
explosive  gas  and  shell  fragments  with  the  chemical  energy  of  the  explosive.5  Two  different  blast  wave  solutions 
have  been  proposed,  one  by  Taylor  in  1949  and  another  in  1961  by  Friedman.  These  blast  wave  solutions  serve  to 
determine  the  flow  properties  after  the  initial  detonation.12,13  In  2007,  Frost  performed  experiments  which  tracked 
the  explosive  clouds  created  by  the  detonation  of  a  heterogeneous  mixture  of  Nitromethane  and  steel  particles.6,7  The 
Frost  study  will  provide  a  basis  to  compare  the  numeric  simulation  to  experimental  data.  The  purpose  of  this  study  is 
to  construct  and  analyze  a  numeric  algorithm  based  on  theoretical  methods  which  can  accurately  model  the 
detonation  and  two-phase  flow  of  a  heterogeneous  high  explosive. 

II.  Technical  Approach 


A.  Velocity  Initialization 

The  inert  metal  particles  experience  an  initial  acceleration  during  the  detonation  of  the  high  explosive  followed 
by  further  acceleration  due  to  the  resulting  blast  wave  flow.  It  is  necessary  to  determine  the  velocity  of  the  particles 
and  gas  from  detonation  which  is  used  as  an  initial  condition  in  the  algorithm.  One  way  to  approximate  this 
initialization  is  to  modify  the  Gurney  equations.  The  Gurney  method  has  been  used  to  estimate  the  velocity  profile 
of  explosive  gasses  and  surrounding  shell  case.  These  equations  are  derived  from  an  energy  balance  between  the 
Gurney  energy  of  the  high  explosive  used  and  the  kinetic  energy  of  the  gas  phase  and  shell  fragments.5  The  Gurney 
equations  require  modification  in  order  to  approximate  a  heterogeneous  mixture.  The  original  Gurney  equations 
assume  a  continuous  mass  of  high  explosive  surrounded  by  one  solid  shell.  These  equations  can  be  modified  to 
include  several  extremely  small  layers  of  explosive  mass  and  metal  particles  which,  when  solved  iteratively  from  the 
core  to  the  outer  explosive  edge  provide  a  velocity  profile.  Many  repeated  and  very  thin  explosive  layers  will 
approximate  a  heterogeneous  mixture.  The  result  is  shown  below  in  Eq.  (1)  where  R0  and  Rt  are  the  outer  and  inner 
radii  of  the  layer  in  question.  This  provides  the  basis  for  the  initial  detonation  velocity  of  the  metal  particles  and  the 
initial  gas  phase  velocity  profile. 


X+2 

5 


3RR+5R 


4 


Rn 


+  RaRi+R 
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(1) 


B.  Particle  Drag  and  Heat  Transfer 

Several  drag  laws  and  heat  transfer  laws  were  tested  during  the  construction  of  the  model.  The  equations  that 
were  found  to  be  the  most  accurate  for  this  type  of  experiment  were  outlined  by  Akhatov  and  are  described  below.8 
Each  spherical  particle  experiences  a  drag  force  given  by  Eq.  (2).  This  force  accounts  for  varying  Reynolds  numbers 
given  by  Eq.  (3)  and  a  piecewise  defined  drag  coefficient  that  is  dependent  on  both  the  Reynolds  number  and  the 
particle  volume  concentration  shown  in  Eq.  (4).  The  heat  transfer  law  shown  in  Eq.  (5)  accounts  for  varying  Nusselt 
numbers,  Nu,  by  Eq.  (6)  where  Re  is  the  Reynolds  number  and  Pr  is  the  Prandtl  number.  Reynolds  numbers,  Nusselt 
numbers,  and  drag  coefficients  are  updated  at  each  time  step  for  each  individual  particle. 
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(5) 

(6) 


C.  Taylor’s  Blast  Wave  Solution 

Two  different  blast  wave  solutions  were  analyzed  during  the  construction  of  this  algorithm.  The  simpler  solution 
was  developed  by  Taylor.  This  similarity  solution  creates  position,  velocity,  pressure,  density,  and  temperature 
profiles  over  time  quantified  primarily  by  the  blast  wave  energy  of  the  specific  explosion  and  derived  from  the 
continuity  of  mass,  momentum,  and  energy  laws.  We  apply  Taylor’s  approximate  form  of  this  similarity  solution. 
Defining  the  blast  wave  energy  as  the  chemical  energy  released  during  detonation  of  the  explosive  allows  solving 
for  the  parameter,  A,  in  Eq.  (7)  which  is  used  in  subsequent  calculations.10  Taylor  calculates  the  blast  radius  as  a 
function  of  the  constant,  A ,  and  time  as  in  Eq.  (8). 12  The  zero  subscript  refers  to  the  ambient  condition. 


f  Ajr 

Eb  =  A1  pa  2+0.185)+- - - 

l  r(r  - 1) 


(0.187) 


(7) 


R(t)  =  (2.5  A)215 12'5  (8) 

To  determine  the  gas  phase  velocity  profile,  the  function  Tj  is  the  dimensionless  radius  of  the  blast  wave’s 
interior  equal  to  the  local  radius  divided  by  the  radius  of  the  outer  blast  wave  at  any  given  time  as  shown  in  Eq.  (9). 
Three  other  functions  of  Tj  are  (/) ,  /,  and  \j/  which  represent  the  non-dimensional  velocity,  pressure,  and  density 
profiles  respectively.  Scaling  these  quantities  requires  multiplying  by  appropriate  parameters  as  defined  in  Eqs.  (10) 
through  (15).  For  simplicity,  /3  is  defined  in  Eq.  (16).  With  these  profiles,  temperature  can  also  be  calculated  by 
using  the  assumption  of  a  calorically  perfect  gas  shown  in  Eq.  (17).9’10’11’12 


Jk. 

II 

(9) 

vg  =  AR -3/V(?7) 

(10) 
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(11) 
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p  =  PvYin) 


(14) 


y(ji) 


[(r +  i)  /(r  -  i)]^3/(x+1) 
[{y^-vp-l)iyfr+mi~r) 


p  =  (iY-\)i{r-\) 


T{r1)  =  A2R-%pJtT1)l[yP,¥m 


(15) 

(16) 
(17) 


The  previous  eleven  equations  define  the  behavior  and  thermodynamic  state  of  the  gaseous  flow  after  the  explosion. 

D.  Friedman’s  Blast  Wave  Solution 

The  second  blast  wave  solution  modeled  in  the  algorithm  was  developed  by  Friedman.  This  solution  produces 
position,  velocity,  pressure,  density,  and  temperature  profiles  of  the  explosive  gasses  and  shock  wave.  The  solution 
was  derived  by  a  spatial  modification  of  a  pressurized,  one-dimensional  shock  tube  and  applied  to  a  much  larger 
explosion.  This  solution  differs  from  Taylor’s  in  that  it  uses  pressure  as  opposed  to  blast  wave  energy  to  quantify  the 
initial  strength  of  the  explosion.  Friedman’s  solution  also  incorporates  the  method  of  characteristics  in  its 
derivation.13 

An  important  quantity  to  be  calculated  is  the  shock  wave  velocity  at  the  instant  of  the  explosion.  This  quantity  is 
computed  by  solving  for  the  thermodynamic  variables  in  a  one-dimensional  shock  tube  with  appropriate  pressure 
and  density  conditions.  Solving  Eq.  (19)  iteratively  provides  the  solution  for  the  initial  Mach  number  of  the 
radiating  shock  wave,  labeled  M.13  For  simplicity,  5  is  defined  in  Eq.  (18)  and  is  a  function  of  the  ratio  of  specific 
heat  of  the  gas  phase.  The  subscript  zero  refers  to  the  ambient  conditions,  four  refers  to  the  conditions  inside  the 
detonated  explosive,  and  i  refers  to  the  initial  time  condition. 

S  =  (y-\)l(y  +  \)  (18) 


{>|  {S+\)/8 

1_— (V.-M,-1)  =(1  +  S)Mi2-S  (19) 

a4  J 

Unlike  a  one-dimensional  shock  tube  where  the  shock  wave  Mach  number  stays  constant  over  time,  the  shock 
wave  from  a  spherical  explosion  will  dampen  significantly.  Friedman  accounts  for  this  slowing  of  the  shock  wave 
with  a  spatial  correction  factor  that  reduces  to  Eq.  (20)  shown  below.13  The  subscript,  s,  refers  to  the  shock  wave. 


2^  =  dM\ 
x„ 


4  M 


■  + 


2  (M2  + 1) 


2/M2  -  y  +  \  m^\2/M2  —/  +  1][M2(/  - 1)  +  2] 


+  ' 


2  M  (zt-l)M2  +2  M2+ 1 


+  - 


M  -l\  2yM  -y  +  \  M(MZ  - 1) 


(20) 


The  shock  wave  position  and  velocity  information  over  time  allows  the  calculation  of  gas  phase  velocity, 
pressure,  density,  and  speed  of  sound  directly  behind  the  shock  using  the  Rankine-Hugoniot  relations  shown  in  Eqs. 
(21),  (22),  (23),  and  (24). 13,15,16  All  of  the  following  thermodynamic  states  remain  constant  in  the  flow  between  the 
shock  wave  and  the  advancing  front  for  the  explosive  gas. 
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(21) 


7  +  1 


P  =  -^—{2yM2  -r  +  l) 

7  +  1 

Po(7  +  1)M2 

P  (7-l)M2+2 


(22) 


(23) 


2  a02(27M2-7  +  l)[M2(7-l)  +  2] 

(7  +  l)2M2 


(24) 


The  border  between  the  advancing  front  of  the  gaseous  detonation  products  and  the  surrounding  air  is  referred  to 
as  the  contact  discontinuity.  It  is  important  to  determine  the  location  of  the  contact  discontinuity  because  at  that 
point  temperature,  density,  and  entropy  experience  an  abrupt  change  even  though  pressure  and  flow  velocity  remain 
constant.  Being  dependent  on  shock  wave  position  and  Mach  number,  the  contact  discontinuity  position  over  time  is 
calculated  with  Eqs.  (25)  and  (26). 13  Eqs.  (27)  and  (28)  define  w  and  Q  which  are  variables  used  in  the  method  of 
characteristics.  The  subscript,  c ,  refers  to  the  contact  discontinuity  while  the  primes  represent  differentiation  with 
respect  to  Mach  number. 


dtq 

dM 


ts)-wts'}/(2Q-w) 


(25) 


xc=xs+w(tc-ts)  (26) 

w  =  vg- a  (27) 

Q  =  a!(y-  l)  +  vg/2  (28) 


Another  important  aspect  of  the  explosion  to  model  is  what  happens  to  the  flow  inside  the  sphere  as  the  high 
pressure  gas  starts  to  travel  outward.  When  this  happens,  an  inward  moving  expansion  wave  is  created  where  gas 
velocity,  speed  of  sound,  pressure,  density,  and  temperature  vary  continuously.  The  algorithm  computes  these 
quantities  in  a  similar  manner  to  a  pressurized  shock  tube  using  Eqs.  (29)  through  (33)  shown  below.14,16,17 

Combined,  these  equations  represent  another  way  to  define  the  behavior  and  flow  of  the  gaseous  detonation 
products.  The  Friedman  blast  wave  solution  is  based  on  a  variety  of  different  assumptions  and  methods  as  compared 
to  the  Taylor  blast  wave  solution  and  thus  produces  significantly  different  results. 
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E.  Energy  Coupling 

The  algorithm  also  incorporates  energy  coupling  into  both  the  Taylor  and  Friedman  blast  wave  solutions.  When 
using  the  Taylor  blast  wave  solution,  the  change  in  particle  kinetic  and  thermal  energy  is  subtracted  from  the  blast 
wave  energy  for  each  particle  per  time  step.  This  results  in  a  general  attenuation  of  the  blast  wave  as  the  particles 
drain  energy.  As  time  progresses  and  the  increased  particle  velocities  and  temperatures  approach  that  of  the  blast 
wave,  less  transfer  occurs.  Thus,  Taylor’s  blast  wave  solution  must  be  re-evaluated,  accounting  for  the  updated  blast 
wave  energy,  at  each  time  step  in  the  algorithm. 

Since  the  Friedman  blast  wave  solution  is  derived  from  the  pressurized  one-dimensional  shock  tube  solution,  it  is 
necessary  to  reduce  the  initial  pressure  condition  in  order  to  model  the  energy  loss  from  the  gaseous  flow.  Energy 
loss  is  determined  at  each  new  time  step  in  the  algorithm  followed  by  an  appropriate  reduction  of  the  initial 
detonation  pressure.  Like  the  Taylor  solution,  the  entire  Friedman  blast  wave  solution  must  be  re-evaluated  at  each 
time  step  in  the  algorithm  to  account  for  energy  transfer  between  the  gas  phase  and  solid  phase. 

F.  Particle  Tracking 

The  gas  phase  and  solid  phase  are  tracked  using  an  Eulerian  and  Lagrangian  approach  respectively.  The  gas 
phase  is  assumed  to  be  a  continuum.  Gas  flow  properties  change  gradually  when  using  the  Taylor  method.  When 
using  the  Friedman  method,  the  flow  properties  change  gradually  as  well  with  the  exceptions  of  the  abrupt  changes 
occurring  at  the  head  of  the  expansion  wave,  tail  of  the  expansion  wave,  and  contact  discontinuity.  Each  particle’s 
position,  velocity,  and  temperature  is  tracked  individually  in  a  Lagrangian  manner  and  updated  at  each  time  step. 

III.  Results  and  Analysis 


A.  Gurney  Velocity  Initialization 

The  results  of  the  Gurney  velocity  initialization  are  shown  in  Fig.  1.  Illustrated  is  increasing  velocity  from  the 
explosive  core  to  the  outer  edge  of  the  explosive.  Gas  near  the  outer  edge  of  the  outer  shell  experiences  less 
backward  pressure  and  thus  achieves  a  higher  initial  velocity.  Particle  velocities  are  slower  than  gas  velocities  at 
each  point  due  to  the  fact  that  each  calculated  layer  of  explosive  has  only  a  fraction  of  surface  area  made  up  of 
particles.  Detonated  gas  escapes  between  the  gaps  in  the  particles  which  creates  the  velocity  difference  shown.  This 
method  produces  an  acceptable  input  to  the  algorithm. 
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Figure  1.  Velocity  as  a  function  of  radial  distance  from  core. 


B.  Taylor  Method  Results 

Using  Taylor’s  blast  wave  solution  in  the  algorithm  produces  several  important  results.  Figure  2  shows  the 
velocity  profile  of  the  gas  phase  during  five  different  times.  Velocity  is  shown  as  a  function  of  radial  distance  from 
the  core  of  the  explosive.  The  velocity  profile  shows  no  abrupt  changes  and  attenuates  as  it  expands  outward.  Figure 
3  illustrates  the  pressure  profile  for  the  gas  phase  for  equivalent  times  in  the  simulation.  Taylor’s  method  predicts  a 
spike  in  pressure  near  the  edge  of  the  propagating  explosive  gas  which  is  produced  by  the  algorithm. 


Distance,  m 

Figure  2.  Gas  velocity  as  a  function  of  distance. 


Distance,  m 

Figure  3.  Gas  pressure  as  a  function  of  distance. 


Figure  4  plots  the  temperature  of  the  gas  phase  as  a  function  of  radial  distance.  This  result  clearly  illustrates  one 
of  the  limitations  associated  with  using  the  Taylor  method,  i.e.,  temperatures  near  the  core  are  far  too  high.  The 
primary  reason  for  this  difficulty  is  that  the  derivation  of  the  temperature  profile  assumes  a  calorically  perfect  gas 
with  a  constant  ratio  of  specific  heats.  This  assumption  is  less  valid  near  the  core  of  the  explosive  where  the  ratio  of 
specific  heats  changes  dramatically.  Near  the  outer  edge  of  the  gas  flow,  a  calorically  perfect  gas  is  more  accurate 
and  yields  more  accurate  temperature  results. 
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Figure  5  displays  gas  density  as  a  function  of  distance  from  the  core  over  various  times.  The  density  profile  is 
similar  to  that  of  the  pressure  profile.  This  physically  means  that  the  bulk  of  the  gaseous  detonation  products 
propagate  outward  in  a  dense  high  pressure  shell.  This  result  shows  another  limitation  for  the  Taylor  blast  wave 
solution  resulting  from  the  assumption  that  the  entire  explosive  has  radiated  outward  starting  from  a  very  small 
source.  The  algorithm  compensates  by  initializing  the  blast  radius  at  the  same  size  as  the  explosive  sphere. 
Nevertheless,  the  assumption  of  a  very  small  source  at  the  initial  time  significantly  affects  the  results  of  the  density 
and  pressure  profiles. 


Distance,  m 

Figure  4.  Gas  temperature  as  a  function  of  distance. 


Distance,  m 

Figure  5.  Gas  density  as  a  function  of  distance. 


Figures  6  and  7  show  the  results  of  the  solid  phase  using  Taylor’s  method.  In  Fig.  6,  particle  velocity  is  slow 
near  the  core  of  the  explosive  and  higher  near  the  outer  edge.  All  particles  slow  over  time  as  they  travel  outward. 
From  Fig.  7,  it  can  be  seen  that  the  particles  heat  rapidly  soon  after  detonation  and  then  gradually  cool.  The 
inaccuracy  in  the  calorically  perfect  gas  model  can  be  seen  from  these  results  as  well.  Several  particles  achieve 
much  greater  temperatures  than  expected  near  the  core. 


Figure  6.  Particle  velocity  as  a  function  of  distance.  Figure  7.  Particle  temperature  as  a  function  of  distance. 


The  validity  of  the  algorithm  can  be  judged  by  comparison  with  the  Frost  data.6  Figures  8  and  9  show  an  overlay 
of  trajectory  plots  created  by  the  algorithm  as  well  as  experimental  data  measured  by  Frost.  Figure  8  graphs  the 
leading  edge  of  the  gas  phase  and  leading  particles  sized  at  463  micrometers  in  diameter.  Figure  9  shows  the  results 
using  925  micrometer  diameter  particles.  All  results  use  nitromethane  as  the  high  explosive. 
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A  crossing  of  the  trajectories  is  shown  on  both  graphs.  Physically,  this  represents  how  the  particles  gradually 
outrun  the  shock  wave.  The  particles  accelerate  during  the  initial  detonation  and  then  accelerate  further  within  the 
resulting  blast  wave  flow.  Over  time,  the  particles  gain  velocity  and  momentum  while  the  shock  wave  starts  to 
dampen  out.  This  effect  is  observed  in  both  the  algorithm  as  well  as  experiment. 


Distance,  m 


Figure  8.  Time  versus  distance  trajectory  plot  including  Frost  data  for  463  micrometer  particles. 


Figure  9.  Time  versus  distance  trajectory  plot  including  Frost  data  for  925  micrometer  particles. 

C.  Friedman  Method  Results 

Using  Friedman’s  blast  wave  solution  in  the  algorithm  produces  distinctly  different  results  than  the  Taylor 
method.  Because  this  method  analyzes  specific  sections  of  flow  separately,  there  are  discontinuous  thermodynamic 
variables.  Velocity  is  shown  in  Fig.  10  to  have  a  distinct  spike  toward  the  outer  shell  and  gradually  slow  near  the 
core.  Pressure  is  maintained  near  the  core  of  the  explosive  rather  than  accumulating  near  the  outer  edge  of  the 
propagating  sphere  displayed  in  Fig.  11  in  a  more  accurate  representation  of  the  flow. 
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Distance,  m  Distance,  m 

Figure  10.  Gas  velocity  as  a  function  of  distance.  Figure  11.  Gas  pressure  as  a  function  of  distance. 


The  temperature  of  the  gas  phase  is  dramatically  different  with  Friedman’s  solution  than  Taylor’s.  The  gas  can 
be  seen  in  Fig.  12  to  be  less  hot  near  the  core  and  cools  much  less  significantly  with  radial  distance.  There  is  no 
asymptotic  increase  in  temperature  near  the  core  as  there  was  with  the  Taylor  solution.  Temperature  is  also  shown  to 
be  higher  in  the  expansion  wave  than  in  the  rest  of  the  flow.  Gas  density  is  seen  to  be  high  near  the  core  in  the 
expansion  wave  and  gradually  decreases  with  radial  distance  as  shown  in  Fig.  13. 


Figure  12.  Gas  temperature  as  a  function  of  distance.  Figure  13.  Density  as  a  function  of  distance. 


The  solid  phase  shows  expected  results.  The  particles  accelerate  very  rapidly  early  in  the  simulation  due  to  the 
initial  detonation  acceleration  and  blast  wave  flow  but  slow  due  to  drag  illustrated  in  Fig.  14.  Figure  15  shows  that 
the  particles  initially  heat  in  the  expansion  wave  but  then  cool  as  they  propagate  outward.  There  are  no  abnormally 
hot  particles  near  the  core  using  the  Friedman  blast  wave  solution. 


10 

American  Institute  of  Aeronautics  and  Astronautics 

092407 


Distance,  m  Distance,  m 

Figure  14.  Particle  velocity  as  a  function  of  distance.  Figure  15.  Particle  temperature  as  a  function  of  distance. 


Figures  16  and  17  show  trajectory  plots  along  with  the  data  from  the  Frost  experiment.6  In  this  case,  the 
Friedman  blast  wave  solution  produces  less  accurate  results  than  the  Taylor  blast  wave  solution.  Here,  the  shock 
wave  propagates  far  too  quickly  at  the  beginning  of  the  simulation  and  thus  propels  the  particles  to  an  unusually 
high  velocity.  The  point  where  the  particles  outrun  the  shock  occurs  too  early  in  time  as  well  as  to  close  to  the  core. 
The  trajectory  plots  expose  an  inaccuracy  in  using  the  Friedman  blast  wave  solution  even  though  the  more  intricate 
nature  of  the  flow  is  captured  compared  to  the  Taylor  solution. 


0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6 

Distance,  m 


Figure  16.  Time  versus  distance  trajectory  plot  including  Frost  data  for  463  micrometer  particles. 
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Figure  17.  Time  versus  distance  trajectory  plot  including  Frost  data  for  925  micrometer  particles. 


IV.  Conclusion 

This  algorithm  has  modeled  a  heterogeneous  high  explosive  by  combining  a  velocity  initialization  scheme, 
appropriate  drag  and  heat  transfer  laws,  two  different  blast  wave  profiles,  kinetic  and  thermal  energy  coupling,  and  a 
Lagrangian-Eulerian  tracking  scheme.  One  of  the  significant  limitations  to  this  model  is  the  assumption  that  the 
detonation  products  act  as  a  calorically  perfect  gas.  To  refine  the  algorithm,  it  is  suggested  to  model  a  thermally 
perfect  gas  which  accounts  for  a  changing  ratio  of  specific  heats.  The  Taylor  blast  wave  solution  is  also  limited  by 
its  assumption  that  the  radiating  gas  products  are  compressed  to  a  small  point  at  the  initial  time.  Still,  Taylor’s  use  of 
a  more  general  blast  wave  energy  quantization  allows  the  method  to  produce  reasonably  accurate  results.  The 
Friedman  blast  wave  solution  does  capture  more  details  of  the  blast  wave  flow;  however  it  produces  less  accurate 
results  displayed  by  its  trajectory  plots.  One  reason  for  the  inaccuracy  while  using  the  Friedman  blast  wave  solution 
is  that  while  the  expansion  wave,  shock  wave,  and  contact  discontinuity  are  modeled  correctly,  flow  properties  at 
later  times  require  additional  modeling.  For  example,  when  the  expansion  wave  reaches  the  core  of  the  explosive,  in 
reality  a  reflected  wave  would  be  generated.  This  is  not  accounted  for  in  Friedman’s  analysis.  Friedman’s  solution 
was  not  based  on  the  extremely  high  pressures  commonly  occurring  in  the  detonation  of  high  explosives.  We  plan  to 
continue  improving  our  implementation  of  the  algorithms  as  time  permits.  Still,  both  solutions  used  in  the  algorithm 
provide  insight  into  the  behavior  of  heterogeneous  high  explosives. 
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Introduction 


i 


•  Purpose 

—  To  develop  and  test  an  accurate  multiphase  blast  explosive 
model 

—  Research  results  of  Gurney  initialization  scheme 

*  Objectives 

—  Incorporate  Gurney  method  initialization 
—  Test  various  drag  laws 
—  Test  various  heat  transfer  laws 
—  Model  Taylor  and  Friedman  Blast  Wave  solutions 
—  Examine  energy  coupling 
—  Compare  to  experiment 
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Background 


i 


*  Taylor  versus  Friedman  blast  wave  solutions 

*  Gurney  method  approximates  particle  velocities  for 
outer  shell  fragments 

*  “Spray”  drag  law,  Carlson  and  Hoglund  rocket 
equations,  Vainshtein  and  Akhatov  explosion 
equations 

*  Frost  and  Zhang  experimental  data  from  2007 
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Technical  Approach 


i 


*  Designed  numerical  algorithm 
—  Initialize  velocity 

—  Determine  particle  acceleration  and  temperature 
change 

—  Alter  blast  wave  solution 
—  Update  particle  position  and  temperature 

*  Combined  different  authors’  work  with  most 
reasonable  assumptions 

*  Compared  with  experimental  data 
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Drag  and  Temperature 
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Taylor’s  Governing  Equations 


•  Pressure 

•  Density 

•  Velocity 

•  Blast  Radius 

•  Temperature 
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Friedman  Method 


i 


*  Modeled  as  a  perturbed  shock  tube 

*  Initial  pressure  assigned  as  Chapman-Jouquet 
pressure 

*  Rankine-Hugoniot  relations  for  flow  behind  shock 
wave 

*  Method  of  characteristics  used  for  contact 
discontinuity 

*  Shock  tube  equations  for  expansion  wave 
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Solid  Phase  Velocity 
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*  Taylor  solution 
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*  Friedman  solution 
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Summary/Conclusions/Recommendations 


*  Provides  theoretical  tool  for  study 

*  Important  in  development  of  heterogeneous  high 
explosives 

*  Recommend  model  as  thermally  perfect  gas,  variable 
specific  heats 

*  Improvements  on  secondary  shock  in  Friedman  model 
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